! written by jxzou at 20200623

module periodic_table
 implicit none
 integer, parameter :: period_nelem = 112
 character(len=2), parameter :: period_elem(0:period_nelem) = (/'Bq',&
   'H ', 'He', 'Li', 'Be', 'B ', 'C ', 'N ', 'O ', 'F ', 'Ne', &
   'Na', 'Mg', 'Al', 'Si', 'P ', 'S ', 'Cl', 'Ar', 'K ', 'Ca', &
   'Sc', 'Ti', 'V ', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', &
   'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', 'Rb', 'Sr', 'Y ', 'Zr', &
   'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', &
   'Sb', 'Te', 'I ', 'Xe', 'Cs', 'Ba', 'La', 'Ce', 'Pr', 'Nd', &
   'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', &
   'Lu', 'Hf', 'Ta', 'W ', 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg', &
   'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn', 'Fr', 'Ra', 'Ac', 'Th', &
   'Pa', 'U ', 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf', 'Es', 'Fm', &
   'Md', 'No', 'Lr', 'Rf', 'Db', 'Sg', 'Bh', 'Hs', 'Mt', 'Ds', &
   'Rg', 'Cn'/)
 ! the radii below are read from Gaussian output file
 real(kind=8), parameter :: vdw_radii(period_nelem) = &
 (/1.4430d0, 1.1810d0, 1.2255d0, 1.3725d0, 2.0415d0, &
   1.9255d0, 1.8300d0, 1.7500d0, 1.6820d0, 1.6215d0, &
   1.4915d0, 1.5105d0, 2.2495d0, 2.1475d0, 2.0735d0, &
   2.0175d0, 1.9735d0, 1.9340d0, 1.9060d0, 1.6995d0, &
   1.6475d0, 1.5875d0, 1.5720d0, 1.5115d0, 1.4805d0, &
   1.4560d0, 1.4360d0, 1.4170d0, 1.7475d0, 1.3815d0, &
   2.1915d0, 2.1400d0, 2.1150d0, 2.1025d0, 2.0945d0, &
   2.0705d0, 2.0570d0, 1.8205d0, 1.6725d0, 1.5620d0, &
   1.5825d0, 1.5260d0, 1.4990d0, 1.4815d0, 1.4645d0, &
   1.4495d0, 1.5740d0, 1.4240d0, 2.2315d0, 2.1960d0, &
   2.2100d0, 2.2350d0, 2.2500d0, 2.2020d0, 2.2585d0, &
   1.8515d0, 1.7610d0, 1.7780d0, 1.8030d0, 1.7875d0, &
   1.7735d0, 1.7600d0, 1.7465d0, 1.6840d0, 1.7255d0, &
   1.7140d0, 1.7045d0, 1.6955d0, 1.6870d0, 1.6775d0, &
   1.8200d0, 1.5705d0, 1.5850d0, 1.5345d0, 1.4770d0, &
   1.5600d0, 1.4200d0, 1.3770d0, 1.6465d0, 1.3525d0, &
   2.1735d0, 2.1485d0, 2.1850d0, 2.3545d0, 2.3750d0, &
   2.3825d0, 2.4500d0, 1.8385d0, 1.7390d0, 1.6980d0, &
   1.7120d0, 1.6975d0, 1.7120d0, 1.7120d0, 1.6905d0, &
   1.6630d0, 1.6695d0, 1.6565d0, 1.6495d0, 1.6430d0, &
   1.6370d0, 1.6240d0, 1.6180d0, 1.7500d0, 1.7500d0, &
   1.7500d0, 1.7500d0, 1.7500d0, 1.7500d0, 1.7500d0, &
   1.7500d0, 1.7500d0/)
contains

 ! map a nuclear charge to an element (e.g. 6->'C')
 ! Note: only 1-112 elements are supported!
 pure function nuc2elem(i) result(s)
  implicit none
  integer, intent(in) :: i
  character(len=2) :: s
 
  s = period_elem(i)
 end function nuc2elem

 pure function elem2vdw_radii(elem) result(radii)
  implicit none
  real(kind=8) :: radii
  character(len=2), intent(in) :: elem

  radii = vdw_radii(elem2nuc(elem))
 end function elem2vdw_radii

 ! map an element to a nuclear charge (e.g. 'C'->6)
 ! Note: only 1-112 elements are supported!
 pure function elem2nuc(s) result(i)
  implicit none
  integer :: i
  character(len=2), intent(in) :: s

  do i = 1, period_nelem, 1
   if(period_elem(i) == s) return
  end do ! for i
 end function elem2nuc

 ! read elements array from a given .gjf array
 subroutine read_elem_from_gjf(gjfname, natom, elem, ghost)
  implicit none
  integer :: i, j, k, nblank, fid
  integer, intent(in) :: natom
  character(len=6) :: str
  character(len=240) :: buf
  character(len=240), intent(in) :: gjfname
  character(len=2), intent(out) :: elem(natom)
  logical, intent(out) :: ghost(natom)

  forall(i = 1:natom) elem(i) = '  '
  ghost = .false.
  open(newunit=fid,file=TRIM(gjfname),status='old',position='rewind')
  nblank = 0

  do while(.true.)
   read(fid,'(A)') buf
   if(LEN_TRIM(buf) == 0) nblank = nblank + 1
   if(nblank == 2) exit
  end do ! for while

  read(fid,'(A)') buf ! skip charge and mult

  do i = 1, natom, 1
   read(fid,*) str
   str = ADJUSTL(str)
   k = LEN_TRIM(str)
   if(str(k-2:k) == '-Bq') then
    elem(i) = str(1:k-3)
    ghost(i) = .true.
   else if(str(2:3)=='(f' .or. str(2:3)=='(F') then
    call upper(str(1:1))
    elem(i) = str(1:1)//' '
   else if(str(3:4)=='(f' .or. str(3:4)=='(F') then
    call upper(str(1:1))
    call lower(str(2:2))
    elem(i) = str(1:2)
   else
    k = IACHAR(str(1:1))
    if((k>64 .and. k<91) .or. (k>96 .and. k<123)) then ! A-Z, a-z
     if(k>96 .and. k<123) str(1:1) = ACHAR(k-32)
     j = IACHAR(str(2:2))
     if(j>64 .and. j<91) str(2:2) = ACHAR(j+32)
     elem(i) = TRIM(str)
    else
     read(str,*) j
     if(j > period_nelem) then
      write(6,'(A)') 'ERROR in subroutine read_elem_from_gjf: j too large.'
      write(6,'(A,I0)') 'j=', j
      stop
     end if
     elem(i) = nuc2elem(j)
    end if
   end if
  end do ! for natom

  close(fid)
 end subroutine read_elem_from_gjf
end module periodic_table

! find the number of atoms in Gaussian .gjf file
subroutine read_natom_from_gjf(gjfname, natom)
 implicit none
 integer :: i, fid, nblank
 integer, intent(out) :: natom
 character(len=240) :: buf
 character(len=240), intent(in) :: gjfname

 nblank = 0
 open(newunit=fid,file=TRIM(gjfname),status='old',position='rewind')
 do while(.true.)
  read(fid,'(A)',iostat=i) buf
  if(i /= 0) exit
  if(LEN_TRIM(buf) == 0) nblank = nblank + 1
  if(nblank == 2) exit
 end do ! for while

 if(i /= 0) then
  write(6,'(A)') 'ERROR in subroutine read_natom_from_gjf: incomplete file '//TRIM(gjfname)
  stop
 end if

 read(fid,'(A)') buf ! skip charge and mult

 natom = 0
 do while(.true.)
  read(fid,'(A)',iostat=i) buf
  if(i /= 0) exit
  if(LEN_TRIM(buf) == 0) exit
  natom = natom + 1
 end do ! for while

 close(fid)
end subroutine read_natom_from_gjf

! read the number of atoms from a given .fch file
subroutine read_natom_from_fch(fchname, natom)
 implicit none
 integer :: i, fid
 integer, intent(out) :: natom
 character(len=240) :: buf
 character(len=240), intent(in) :: fchname

 natom = 0
 open(newunit=fid,file=TRIM(fchname),status='old',position='rewind')

 do while(.true.)
  read(fid,'(A)',iostat=i) buf
  if(buf(1:14) == 'Atomic numbers') exit
 end do ! for while

 BACKSPACE(fid)
 read(fid,'(A49,2X,I10)') buf, natom

 close(fid)
end subroutine read_natom_from_fch

! read the number of atoms from a .xyz file
subroutine read_natom_from_xyz(xyzname, natom)
 implicit none
 integer :: i, fid
 integer, intent(out) :: natom
 character(len=240), intent(in) :: xyzname

 natom = 0
 open(newunit=fid,file=TRIM(xyzname),status='old',position='rewind')
 read(fid,*,iostat=i) natom
 close(fid)

 if(i /= 0) then
  write(6,'(A)') 'ERROR in subroutine read_natom_from_xyz: failed to read natom&
                 & from file '//TRIM(xyzname)
  stop
 end if
end subroutine read_natom_from_xyz

! read the number of atoms from a .pdb file
! Note: if there exists >1 frames, only the first frame will be dectected
subroutine read_natom_from_pdb(pdbname, natom)
 implicit none
 integer :: i, fid
 integer, intent(out) :: natom
 character(len=13) :: buf
 character(len=240), intent(in) :: pdbname

 natom = 0
 open(newunit=fid,file=TRIM(pdbname),status='old',position='rewind')
 do while(.true.)
  read(fid,'(A)',iostat=i) buf
  if(i /= 0) exit
  if(buf(1:3) == 'END') exit
 end do ! for while

 if(i /= 0) then
  write(6,'(A)') 'ERROR in subroutine read_natom_from_pdb: failed to read&
                   & natom from file '//TRIM(pdbname)
  close(fid)
  stop
 end if

 BACKSPACE(fid)
 BACKSPACE(fid)
 read(fid,'(A)') buf
 if(buf(1:3) == 'TER') then
  BACKSPACE(fid)
  BACKSPACE(fid)
  read(fid,'(A)') buf
 end if
 close(fid)

 i = INDEX(buf, ' ')
 read(buf(i+1:),*) natom
end subroutine read_natom_from_pdb

! read the number of atoms from a (Open)Molcas output file
subroutine read_natom_from_molcas_out(outname, natom)
 implicit none
 integer :: i, fid
 integer, intent(out) :: natom
 character(len=240) :: buf
 character(len=240), intent(in) :: outname

 natom = 0
 open(newunit=fid,file=TRIM(outname),status='old',position='rewind')

 do while(.true.)
  read(fid,'(A)',iostat=i) buf
  if(i /= 0) exit
  if(buf(1:21) == '++    Molecular struc') exit
 end do ! for while

 if(i /= 0) then
  write(6,'(A)') "ERROR in subroutine read_natom_from_molcas_out: keywords&
                   & '++    Molecular struc' not found"
  write(6,'(A)') 'in file '//TRIM(outname)
  close(fid)
  stop
 end if

 do while(.true.)
  read(fid,'(A)') buf
  if(buf(6:11) == 'Center') exit 
 end do ! for while

 do while(.true.)
  read(fid,'(A)') buf
  if(LEN_TRIM(buf) == 0) exit
  natom = natom + 1 
 end do ! for while

 close(fid)
end subroutine read_natom_from_molcas_out

! read the number of atoms from a Molpro output file
subroutine read_natom_from_molpro_out(outname, natom)
 implicit none
 integer :: i, fid
 integer, intent(out) :: natom
 character(len=240) :: buf
 character(len=240), intent(in) :: outname

 natom = 0
 open(newunit=fid,file=TRIM(outname),status='old',position='rewind')

 do while(.true.)
  read(fid,'(A)',iostat=i) buf
  if(i /= 0) exit
  if(buf(2:12) == 'ATOMIC COOR') exit
 end do ! for while

 if(i /= 0) then
  write(6,'(A)') "ERROR in subroutine read_natom_from_molpro_out: no '&
                   &ATOMIC COOR' found in file "//TRIM(outname)
  close(fid)
  stop
 end if

 do i = 1, 3
  read(fid,'(A)') buf
 end do

 do while(.true.)
  read(fid,'(A)') buf
  if(LEN_TRIM(buf) == 0) exit
  natom = natom + 1
 end do ! for while

 close(fid)
end subroutine read_natom_from_molpro_out

! read the number of atoms from ORCA .engrad file
subroutine read_natom_from_engrad(fname, natom)
 implicit none
 integer :: fid
 integer, intent(out) :: natom
 character(len=240) :: buf
 character(len=240), intent(in) :: fname

 natom = 0
 open(newunit=fid,file=TRIM(fname),status='old',position='rewind')
 read(fid,'(A)') buf
 read(fid,'(A)') buf

 if(buf(1:16) /= '# Number of atom') then
  write(6,'(/,A)') "ERROR in subroutine read_natom_from_engrad: '# Number of at&
                   &om' expected, but"
  write(6,'(A)') "got '"//TRIM(buf)//"'"
  close(fid)
  stop
 end if

 read(fid,'(A)') buf
 read(fid,*) natom
 close(fid)
end subroutine read_natom_from_engrad

subroutine read_coor_from_gjf_or_xyz(fname, natom, coor)
 implicit none
 integer :: i, charge, mult
 integer, intent(in) :: natom
 integer, allocatable :: nuc(:)
 real(kind=8), intent(out) :: coor(3,natom)
 character(len=2), allocatable :: elem(:)
 character(len=240), intent(in) :: fname

 i = INDEX(fname, '.xyz', back=.true.)
 if(i > 0) then
  allocate(elem(natom))
  call read_elem_and_coor_from_xyz(fname, natom, elem, coor)
  deallocate(elem)
 else ! not .xyz, assume it be .gjf
  allocate(elem(natom), nuc(natom))
  call read_elem_and_coor_from_gjf(fname, natom, elem, nuc, coor, charge, mult)
  deallocate(elem, nuc)
 end if
end subroutine read_coor_from_gjf_or_xyz

! read 3 arrays elem, nuc, coor, and the total charge as well as multiplicity
! from a given .gjf file
subroutine read_elem_and_coor_from_gjf(gjfname, natom, elem, nuc, coor, charge, mult)
 use phys_cons, only: Bohr_const
 use periodic_table, only: elem2nuc
 implicit none
 integer :: i, j, k, fid, nblank, ne
 integer, intent(in) :: natom
 integer, intent(out) :: nuc(natom), charge, mult
 real(kind=8), intent(out) :: coor(3,natom)
 character(len=2), intent(out) :: elem(natom)
 character(len=240) :: buf
 character(len=240), intent(in) :: gjfname
 logical :: bohr

 charge = 0; mult = 1; nuc = 0; coor = 0d0; bohr = .false.; nblank = 0

 open(newunit=fid,file=TRIM(gjfname),status='old',position='rewind')
 do while(.true.)
  read(fid,'(A)',iostat=i) buf
  if(i /= 0) exit
  if(LEN_TRIM(buf) == 0) nblank = nblank + 1
  if(nblank == 2) exit

  if(buf(1:1) =='#') then
   if(INDEX(buf,'units=au') /= 0) bohr = .true.
  end if
 end do ! for while

 if(i /= 0) then
  write(6,'(A)') 'ERROR in subroutine read_elem_and_coor_from_gjf: incomplete&
                & file '//TRIM(gjfname)
  stop
 end if

 read(fid,*,iostat=k) charge, mult
 if(k /= 0) then
  write(6,'(A)') 'ERROR in subroutine read_elem_and_coor_from_gjf: failed to&
                & read charge and mult.'
  write(6,'(A)') 'There exists syntax error in file '//TRIM(gjfname)
  stop
 end if

 do i = 1, natom, 1
  read(fid,'(A)') buf
  j = INDEX(buf,'('); k = INDEX(buf,')') ! in case for the fragment guess wfn
  if(j*k /= 0) buf(j:k) = ' '

  read(buf,*,iostat=k) elem(i), coor(1:3,i)
  if(k /= 0) then
   write(6,'(A)') 'ERROR in subroutine read_elem_and_coor_from_gjf: only 4-column&
                 & format is supported.'
   close(fid)
   stop
  end if
 end do ! for i

 close(fid)
 if(bohr) coor = coor*Bohr_const ! convert Bohr to Angstrom

 ! standardize a set of elements, e.g. he -> He
 call standardize_elem(natom, elem)

 ! convert atomic numbers to atomic symbols if the user uses atomic numbers
 call possible_nuc2elem(natom, elem)

 ! convert element symbols to atomic order
 forall(i = 1:natom) nuc(i) = elem2nuc(elem(i))

 ne = SUM(nuc) - charge
 if(MOD(ne,2) /= MOD(mult-1,2)) then
  write(6,'(/,A)') 'ERROR in subroutine read_elem_and_coor_from_gjf:'
  write(6,'(2(A,I0),A)') 'The combination of multiplicity ', mult, ' and ', ne,&
                         &' electrons is impossible.'
  stop
 end if
end subroutine read_elem_and_coor_from_gjf

! convert atomic numbers to atomic symbols if the user uses atomic numbers
subroutine possible_nuc2elem(natom, elem)
 use periodic_table, only: period_elem
 implicit none
 integer :: i, j
 integer, intent(in) :: natom
 character(len=2), intent(inout) :: elem(natom)

 do i = 1, natom, 1
  j = IACHAR(elem(i)(1:1))
  if(j>47 .and. j<58) then
   read(elem(i),*) j
   elem(i) = period_elem(j)
  end if
 end do ! for i
end subroutine possible_nuc2elem

! read charge, spin multiplicities and atom2frag from a given .gjf file
subroutine read_frag_guess_from_gjf(gjfname, natom, atom2frag, nfrag, frag_char_mult)
 implicit none
 integer :: i, j, k, nblank, charge, mult, fid
 integer, intent(in) :: natom, nfrag
 integer, intent(out) :: atom2frag(natom), frag_char_mult(2,nfrag)
 character(len=240) :: buf
 character(len=240), intent(in) :: gjfname

 atom2frag = 0; frag_char_mult = 0
 open(newunit=fid,file=TRIM(gjfname),status='old',position='rewind')

 nblank = 0
 do while(.true.)
  read(fid,'(A)') buf
  if(LEN_TRIM(buf) == 0) nblank = nblank + 1
  if(nblank == 2) exit
 end do ! for while

 read(fid,*,iostat=k) charge, mult, ((frag_char_mult(j,i),j=1,2),i=1,nfrag)
 if(k /= 0) then
  write(6,'(A)') 'ERROR in subroutine read_frag_guess_from_gjf: failed to read&
                & charges and spin'
  write(6,'(A)') 'multiplicities of fragments. Please check syntax in file '//&
                 TRIM(gjfname)
  close(fid)
  stop
 end if

 do i = 1, natom, 1
  read(fid,'(A)') buf
  j = INDEX(buf,'='); k = INDEX(buf,')')

  if(j*k == 0) then
   write(6,'(A)') 'ERROR in subroutine read_frag_guess_from_gjf: failed to read&
                  & atom2frag.'
   write(6,'(A)') 'Problematic line: '//TRIM(buf)
   write(6,'(A)') 'Please check syntax in file '//TRIM(gjfname)
   close(fid)
   stop
  end if

  read(buf(j+1:k-1),*) atom2frag(i)
 end do ! for i

 close(fid)
end subroutine read_frag_guess_from_gjf

! read nuclear charge number from a given .fch file
subroutine read_nuc_from_fch(fchname, natom, nuc)
 implicit none
 integer :: i, fid
 integer, intent(in) :: natom
 integer, intent(out) :: nuc(natom)
 character(len=240) :: buf
 character(len=240), intent(in) :: fchname

 nuc = 0
 open(newunit=fid,file=TRIM(fchname),status='old',position='rewind')

 do while(.true.)
  read(fid,'(A)') buf
  if(buf(1:14) == 'Atomic numbers') exit
 end do ! for while

 read(fid,'(6(1X,I11))') (nuc(i),i=1,natom)
 close(fid)
end subroutine read_nuc_from_fch

! read the Cartesian coordinates from a given .fch file
subroutine read_coor_from_fch(fchname, natom, coor)
 use phys_cons, only: Bohr_const
 implicit none
 integer :: i, fid
 integer, intent(in) :: natom
 real(kind=8), intent(out) :: coor(3,natom)
 real(kind=8), allocatable :: coor0(:)
 character(len=240) :: buf
 character(len=240), intent(in) :: fchname

 open(newunit=fid,file=TRIM(fchname),status='old',position='rewind')

 ! find and read coordinates
 do while(.true.)
  read(fid,'(A)',iostat=i) buf
  if(i /= 0) exit
  if(buf(1:12) == 'Current cart') exit
 end do

 if(i /= 0) then
  write(6,'(A)') "ERROR in subroutine read_coor_from_fch: no 'Current &
                   & cart' found in file "//TRIM(fchname)
  close(fid)
  stop
 end if

 allocate(coor0(3*natom), source=0d0)
 read(fid,'(5(1X,ES15.8))') (coor0(i),i=1,3*natom)
 close(fid)

 coor0 = coor0*Bohr_const ! convert Bohr to Angstrom
 coor = RESHAPE(coor0,(/3,natom/))
 deallocate(coor0)
end subroutine read_coor_from_fch

! read 3 arrays elem, nuc, coor, and the total charge as well as multiplicity
! from a given .fch file
subroutine read_elem_and_coor_from_fch(fchname, natom, elem, nuc, coor, charge, mult)
 use phys_cons, only: Bohr_const
 use periodic_table, only: nuc2elem
 implicit none
 integer :: i, fid
 integer, intent(in) :: natom
!f2py intent(in) :: natom
 integer, intent(out) :: charge, mult, nuc(natom)
 real(kind=8), intent(out) :: coor(3,natom)
 real(kind=8), allocatable :: coor0(:)
 character(len=2), intent(out) :: elem(natom)
!f2py intent(out) :: elem, nuc, coor, charge, mult
!f2py depend(natom) :: elem, nuc, coor
 character(len=240) :: buf
 character(len=240), intent(in) :: fchname
!f2py intent(in) :: fchname

 open(newunit=fid,file=TRIM(fchname),status='old',position='rewind')

 ! find charge and mult
 do while(.true.)
  read(fid,'(A)') buf
  if(buf(1:6) == 'Charge') exit
 end do
 BACKSPACE(fid)
 read(fid,'(A49,2X,I10)') buf, charge
 read(fid,'(A49,2X,I10)') buf, mult

 ! find atomic numbers/nuclear charges
 do while(.true.)
  read(fid,'(A)') buf
  if(buf(1:14) == 'Atomic numbers') exit
 end do ! for i
 read(fid,'(6(1X,I11))') (nuc(i),i=1,natom)

 ! find and read coordinates
 do while(.true.)
  read(fid,'(A)') buf
  if(buf(1:12) == 'Current cart') exit
 end do
 allocate(coor0(3*natom), source=0d0)
 read(fid,'(5(1X,ES15.8))') (coor0(i),i=1,3*natom)
 coor = RESHAPE(coor0,[3,natom])

 deallocate(coor0)
 close(fid)

 coor = coor*Bohr_const ! convert Bohr to Angstrom
 forall(i=1:natom) elem(i) = nuc2elem(nuc(i))
end subroutine read_elem_and_coor_from_fch

! calculate the nuclear dipole moment from
subroutine get_nuc_dipole(natom, nuc, coor, n_dipole)
 use phys_cons, only: Bohr_const
 implicit none
 integer :: i
 integer, intent(in) :: natom
!f2py intent(in) :: natom
 integer, intent(in) :: nuc(natom)
 real(kind=8), allocatable :: rnuc(:)
 real(kind=8), intent(in) :: coor(3,natom)
!f2py intent(in) :: nuc, coor
!f2py depend(natom) :: nuc, coor
 real(kind=8), intent(out) :: n_dipole(3) ! x,y,z 3-components
!f2py intent(out) :: n_dipole

 allocate(rnuc(natom))
 forall(i = 1:natom) rnuc(i) = DBLE(nuc(i))

 do i = 1, 3
  n_dipole(i) = DOT_PRODUCT(rnuc, coor(i,:))
 end do ! for i
 deallocate(rnuc)

 ! input coor are in Angstrom, convert n_dipole into a.u.
 n_dipole = n_dipole/Bohr_const
end subroutine get_nuc_dipole

! read Cartesian gradient from a given file
subroutine read_grad_from_output(prog_name, outname, natom, grad)
 implicit none
 integer, intent(in) :: natom
 real(kind=8) :: e
 real(kind=8), intent(out) :: grad(3*natom)
 character(len=10), intent(in) :: prog_name
 character(len=240), intent(in) :: outname

 select case(prog_name)
 case('bdf')
  call read_grad_from_bdf_out(outname, natom, grad)
 case('cfour')
  call read_grad_from_cfour_out(outname, natom, grad)
 case('dalton')
  call read_grad_from_dalton_out(outname, natom, grad)
 case('gamess')
  call read_grad_from_gms_dat(outname, natom, grad)
 case('gaussian')
  call read_grad_from_gau_log(outname, natom, grad)
 case('molpro')
  call read_grad_from_molpro_out(outname, natom, grad)
 case('orca')
  call read_grad_from_engrad(outname, natom, e, grad)
 case('psi4')
  call read_grad_from_psi4_out(outname, natom, grad)
 case('pyscf')
  call read_grad_from_pyscf_out(outname, natom, grad)
 case('qchem')
  call read_grad_from_qchem131(natom, grad, .true.)
 case('openmolcas')
  call read_grad_from_molcas_out(outname, natom, grad)
 case default
  write(6,'(/,A)') 'ERROR in subroutine read_grad_from_output: program cannot b&
                   &e identified.'
  write(6,'(A)') 'prog_name='//TRIM(prog_name)
  stop
 end select

 write(6,'(/,A)') 'Cartesian gradients (HARTREE/BOHR):'
 write(6,'(5(1X,ES15.8))') grad
end subroutine read_grad_from_output

! read Cartesian gradient from a given PySCF output file
subroutine read_grad_from_pyscf_out(outname, natom, grad)
 implicit none
 integer :: i, k, fid
 integer, intent(in) :: natom
 real(kind=8), intent(out) :: grad(3*natom)
 character(len=3) :: elem
 character(len=240) :: buf
 character(len=240), intent(in) :: outname

 grad = 0d0
 open(newunit=fid,file=TRIM(outname),status='old',position='append')
 do while(.true.)
  BACKSPACE(fid)
  BACKSPACE(fid)
  read(fid,'(A)') buf
  if(INDEX(buf,'gradients') /= 0) exit
 end do ! for while

 read(fid,'(A)') buf

 do i = 1, natom, 1
  read(fid,*) k, elem, grad(3*i-2:3*i)
 end do ! for i

 close(fid)
end subroutine read_grad_from_pyscf_out

! read Cartesian gradient from a given .fch file
subroutine read_grad_from_fch(fchname, natom, grad)
 implicit none
 integer :: i, fid
 integer, intent(in) :: natom
 real(kind=8), intent(out) :: grad(3*natom)
 character(len=240) :: buf
 character(len=240), intent(in) :: fchname

 grad = 0d0
 open(newunit=fid,file=TRIM(fchname),status='old',position='rewind')
 do while(.true.)
  read(fid,'(A)') buf
  if(buf(1:26) == 'Opt point       1 Gradient') exit
 end do ! for while

 read(fid,'(5(1X,ES15.8))') (grad(i),i=1,3*natom)
 close(fid)
end subroutine read_grad_from_fch

! read Cartesian gradient from a given Gaussian .log file
subroutine read_grad_from_gau_log(logname, natom, grad)
 implicit none
 integer :: i, k1, k2, fid
 integer, intent(in) :: natom
 real(kind=8), intent(out) :: grad(3*natom)
 character(len=240) :: buf
 character(len=240), intent(in) :: logname

 grad = 0d0
 open(newunit=fid,file=TRIM(logname),status='old',position='rewind')
 do while(.true.)
  read(fid,'(A)') buf
  if(buf(38:52) == 'Forces (Hartree') exit
 end do ! for while

 read(fid,'(A)') buf
 read(fid,'(A)') buf

 do i = 1, natom, 1
  read(fid,*) k1, k2, grad(3*i-2:3*i)
 end do ! for i

 close(fid)

 grad = -grad !!! VIP
end subroutine read_grad_from_gau_log

! read Cartesian gradient from a given GAMESS .gms file
subroutine read_grad_from_gms_gms(outname, natom, grad)
 implicit none
 integer :: i, k, fid
 integer, intent(in) :: natom
 real(kind=8), intent(out) :: grad(3*natom)
 character(len=2) :: elem = ' '
 character(len=240) :: buf
 character(len=240), intent(in) :: outname

 grad = 0d0

 open(newunit=fid,file=TRIM(outname),status='old',position='append')
 do while(.true.)
  BACKSPACE(fid)
  BACKSPACE(fid)
  read(fid,'(A)') buf
  if(buf(26:47) == 'GRADIENT OF THE ENERGY') exit
 end do ! for while

 do i = 1, 3
  read(fid,'(A)') buf
 end do ! for i

 do i = 1, natom, 1
  read(fid,*) k, elem, grad(3*i-2:3*i)
 end do ! for i

 close(fid)
end subroutine read_grad_from_gms_gms

! read Cartesian gradient from a given GAMESS .dat file
subroutine read_grad_from_gms_dat(datname, natom, grad)
 implicit none
 integer :: i, fid
 integer, intent(in) :: natom
 real(kind=4) :: r
 real(kind=8), intent(out) :: grad(3*natom)
 character(len=2) :: elem
 character(len=240) :: buf
 character(len=240), intent(in) :: datname

 grad = 0d0
 open(newunit=fid,file=TRIM(datname),status='old',position='append')
 do while(.true.)
  BACKSPACE(fid)
  BACKSPACE(fid)
  read(fid,'(A)') buf
  if(buf(2:6) == '$GRAD') exit
 end do ! for while

 read(fid,'(A)') buf
 do i = 1, natom, 1
  read(fid,*) elem, r, grad(3*i-2:3*i)
 end do ! for i

 close(fid)
end subroutine read_grad_from_gms_dat

! read Cartesian gradient from a given (Open)Molcas .out file
subroutine read_grad_from_molcas_out(outname, natom, grad)
 implicit none
 integer :: i, fid
 integer, intent(in) :: natom
 real(kind=8), intent(out) :: grad(3*natom)
 character(len=8) :: str = ' '
 character(len=240) :: buf
 character(len=240), intent(in) :: outname

 grad = 0d0
 open(newunit=fid,file=TRIM(outname),status='old',position='append')

 do while(.true.)
  BACKSPACE(fid)
  BACKSPACE(fid)
  read(fid,'(A)') buf
  if(buf(17:35) == 'Molecular gradients') exit
 end do ! for while

 do i = 1, 7, 1
  read(fid,'(A)') buf
 end do ! for i

 do i = 1, natom, 1
  read(fid,*) str, grad(3*i-2:3*i)
 end do ! for i

 close(fid)
end subroutine read_grad_from_molcas_out

! read Cartesian gradient from a given ORCA .out file
subroutine read_grad_from_orca_out(outname, natom, grad)
 implicit none
 integer :: i, k, fid
 integer, intent(in) :: natom
 real(kind=8), intent(out) :: grad(3*natom)
 character(len=1) :: str = ' '
 character(len=2) :: elem = ' '
 character(len=240) :: buf
 character(len=240), intent(in) :: outname

 grad = 0d0
 open(newunit=fid,file=TRIM(outname),status='old',position='append')

 do while(.true.)
  BACKSPACE(fid)
  BACKSPACE(fid)
  read(fid,'(A)') buf
  if(buf(1:14) == 'CARTESIAN GRAD') exit
 end do ! for while

 read(fid,'(A)') buf
 read(fid,'(A)') buf

 do i = 1, natom, 1
  read(fid,*) k, elem, str, grad(3*i-2:3*i)
 end do ! for i

 close(fid)
end subroutine read_grad_from_orca_out

! read Cartesian gradients from a given ORCA .engrad file
subroutine read_grad_from_engrad(engrad, natom, e, grad)
 implicit none
 integer :: i, fid
 integer, intent(in) :: natom
 real(kind=8), intent(out) :: e, grad(3*natom)
 character(len=240) :: buf
 character(len=240), intent(in) :: engrad

 e = 0d0; grad = 0d0
 open(newunit=fid,file=TRIM(engrad),status='old',position='rewind')

 do while(.true.)
  read(fid,'(A)',iostat=i) buf
  if(i /= 0) exit
  if(buf(1:17) == '# The current tot') exit
 end do ! for while

 if(i /= 0) then
  write(6,'(/,A)') "ERROR in subroutine read_grad_from_engrad: no '# The curren&
                   &t tot' found in"
  write(6,'(A)') 'file '//TRIM(engrad)
  close(fid)
  stop
 end if

 read(fid,'(A)') buf
 read(fid,*) e

 do i = 1, 3
  read(fid,'(A)') buf
 end do
 read(fid,*) grad

 close(fid)
end subroutine read_grad_from_engrad

! read CASSCF or CASPT2 Cartesian gradients from a given Molpro output file
subroutine read_grad_from_molpro_out(outname, natom, grad)
 implicit none
 integer :: i, k, fid
 integer, intent(in) :: natom
 real(kind=8), intent(out) :: grad(3*natom)
 character(len=240) :: buf
 character(len=240), intent(in) :: outname

 grad = 0d0
 open(newunit=fid,file=TRIM(outname),status='old',position='append')

 do while(.true.)
  BACKSPACE(fid)
  BACKSPACE(fid)
  read(fid,'(A)') buf
  if(INDEX(buf,'GRADIENT FOR S') > 0) exit
 end do ! for while

 do i = 1, 3
  read(fid,'(A)') buf
 end do ! for i

 do i = 1, natom, 1
  read(fid,*) k, grad(3*i-2:3*i)
 end do ! for i

 close(fid)
end subroutine read_grad_from_molpro_out

! read Cartesian gradient from a given BDF .out file
subroutine read_grad_from_bdf_out(outname, natom, grad)
 implicit none
 integer :: i, fid
 integer, intent(in) :: natom
 real(kind=8), intent(out) :: grad(3*natom)
 character(len=10) :: str
 character(len=240) :: buf
 character(len=240), intent(in) :: outname

 grad = 0d0
 open(newunit=fid,file=TRIM(outname),status='old',position='append')
 do while(.true.)
  BACKSPACE(fid)
  BACKSPACE(fid)
  read(fid,'(A)') buf
  if(INDEX(buf,'Molecular gradient - Mol') /= 0) exit
 end do ! for while

 read(fid,'(A)') buf
 read(fid,'(A)') buf

 do i = 1, natom, 1
  read(fid,*) str, grad(3*i-2:3*i)
 end do ! for i

 close(fid)
end subroutine read_grad_from_bdf_out

subroutine read_grad_from_psi4_out(outname, natom, grad)
 implicit none
 integer :: i, fid
 integer, intent(in) :: natom
 character(len=10) :: str
 character(len=240) :: buf
 character(len=240), intent(in) :: outname
 real(kind=8), intent(out) :: grad(3*natom)

 grad = 0d0
 open(newunit=fid,file=TRIM(outname),status='old',position='append')

 do while(.true.)
  BACKSPACE(fid,iostat=i)
  if(i /= 0) exit
  BACKSPACE(fid,iostat=i)
  if(i /= 0) exit
  read(fid,'(A)') buf
  if(buf(4:13)=='Total Grad' .or. buf(4:13)=='Total grad') exit
 end do ! for while

 if(i /= 0) then
  write(6,'(/,A)') "ERROR in subroutine read_grad_from_psi4_out: no 'Total Grad&
                   &' found in"
  write(6,'(A)') 'file '//TRIM(outname)
  close(fid)
  stop
 end if

 read(fid,'(A)') buf
 read(fid,'(A)') buf
 do i = 1, natom, 1
  read(fid,*) str, grad(3*i-2:3*i)
 end do ! for i
 close(fid)
end subroutine read_grad_from_psi4_out

subroutine read_grad_from_qchem131(natom, grad, deleted)
 implicit none
 integer :: i, fid
 integer, intent(in) :: natom
 real(kind=8), intent(out) :: grad(3*natom)
 logical, intent(in) :: deleted

 grad = 0d0
 open(newunit=fid,file='131.0',status='old',access='stream')
 read(fid,iostat=i) grad

 if(i /= 0) then
  write(6,'(/,A)') 'ERROR in subroutine read_grad_from_qchem_out: failed to rea&
                   &d nuclear gradients'
  write(6,'(A)') 'from file 131.0'
  close(fid)
  stop
 else ! read successfully
  if(deleted) then
   close(fid,status='delete')
  else
   close(fid)
  end if
 end if
end subroutine read_grad_from_qchem131

subroutine read_grad_from_cfour_out(outname, natom, grad)
 implicit none
 integer :: i, fid
 integer, intent(in) :: natom
 character(len=240) :: buf
 character(len=240), intent(in) :: outname
 real(kind=8), intent(out) :: grad(3*natom)

 grad = 0d0
 open(newunit=fid,file=TRIM(outname),status='old',position='append')

 do while(.true.)
  BACKSPACE(fid,iostat=i)
  if(i /= 0) exit
  BACKSPACE(fid,iostat=i)
  if(i /= 0) exit
  read(fid,'(A)') buf
  if(buf(3:10) == 'gradient') exit
 end do ! for while

 if(i /= 0) then
  write(6,'(/,A)') "ERROR in subroutine read_grad_from_cfour_out: keyword 'grad&
                   &ient' not found"
  write(6,'(A)') 'in file '//TRIM(outname)
  close(fid)
  stop
 end if

 read(fid,*) grad
 close(fid)
end subroutine read_grad_from_cfour_out

subroutine read_grad_from_dalton_out(outname, natom, grad)
 implicit none
 integer :: i, fid
 integer, intent(in) :: natom
 character(len=3) :: str
 character(len=240) :: buf
 character(len=240), intent(in) :: outname
 real(kind=8), intent(out) :: grad(3*natom)

 grad = 0d0
 open(newunit=fid,file=TRIM(outname),status='old',position='append')

 do while(.true.)
  BACKSPACE(fid,iostat=i)
  if(i /= 0) exit
  BACKSPACE(fid,iostat=i)
  if(i /= 0) exit
  read(fid,'(A)') buf
  if(buf(30:43)=='Molecular grad' .or. buf(33:46)=='Molecular grad') exit
 end do ! for while

 if(i /= 0) then
  write(6,'(/,A)') "ERROR in subroutine read_grad_from_dalton_out: keyword 'Mol&
                   &ecular grad' not"
  write(6,'(A)') 'found in file '//TRIM(outname)
  close(fid)
  stop
 end if

 read(fid,'(A)') buf
 read(fid,'(A)') buf

 do i = 1, natom, 1
  read(fid,*) str, grad(3*i-2:3*i)
 end do ! for i

 close(fid)
end subroutine read_grad_from_dalton_out

! read Cartesian xyz coordinates from a .xyz file
! Note: 1) return array coor(3,natom) are in unit Angstrom
!       2) if 'bohr' key is found in the 2nd line of the xyz file,
!          coor will be multiplied by Bohr_const
subroutine read_elem_and_coor_from_xyz(xyzname, natom, elem, coor)
 use phys_cons, only: Bohr_const
 implicit none
 integer :: i, k, fid
 integer, intent(in) :: natom
 real(kind=8), intent(out) :: coor(3,natom)
 character(len=2), intent(out) :: elem(natom)
 character(len=240) :: buf
 character(len=240), intent(in) :: xyzname
 logical :: bohr

 elem = '  '; coor = 0d0
 open(newunit=fid,file=TRIM(xyzname),status='old',position='rewind')
 read(fid,'(A)') buf
 read(fid,'(A)') buf

 bohr = .false.
 call lower(buf)
 if(INDEX(buf,'bohr') > 0) then
  if(INDEX(buf,'angstrom') > 0) then
   write(6,'(A)') "ERROR in subroutine read_elem_and_coor_from_xyz: it's confus&
                  &ing because both"
   write(6,'(A)') "'bohr' and 'angstrom' are detected in the 2nd line of file "&
                  //TRIM(xyzname)
   close(fid)
   stop
  else
   bohr = .true.
  end if
 end if

 do i = 1, natom, 1
  read(fid,*,iostat=k) elem(i), coor(1:3,i)
  if(k /= 0) then
   write(6,'(A)') 'ERROR in subroutine read_elem_and_coor_from_xyz: insufficien&
                  &t number of atoms'
   write(6,'(A)') 'in file '//TRIM(xyzname)
   write(6,'(2(A,I0))') 'Input natom=', natom, ', but broken at i=', i
   close(fid)
   stop
  end if
 end do ! for i

 close(fid)
 if(bohr) coor = coor*Bohr_const ! convert Bohr to Angstrom
end subroutine read_elem_and_coor_from_xyz

! read the number of frames from xyz file
subroutine read_nframe_from_xyz(xyzname, nframe)
 implicit none
 integer :: i, fid, natom
 integer, intent(out) :: nframe
 character(len=240) :: buf
 character*240, intent(in) :: xyzname

 nframe = 0
 open(newunit=fid,file=TRIM(xyzname),status='old',position='rewind')

 do while(.true.)
  read(fid,*,iostat=i) natom
  if(i /= 0) exit
  do i = 1, natom+1, 1
   read(fid,'(A)') buf
  end do ! for i
  nframe = nframe + 1
 end do ! for while

 close(fid)
end subroutine read_nframe_from_xyz

! read the i-th frame from a given .xyz file
subroutine read_iframe_from_xyz(xyzname, ith, natom, elem, coor)
 implicit none
 integer :: i, k, fid, nframe
 integer, intent(in) :: ith, natom
 real(kind=8), dimension(3,natom), intent(out) :: coor
 character(len=240) :: buf
 character*2, dimension(natom), intent(out) :: elem
 character*240, intent(in) :: xyzname

 nframe = 0; coor = 0d0; elem = ' '
 open(newunit=fid,file=TRIM(xyzname),status='old',position='rewind')

 if(ith > 1) then
  do while(.true.)
   read(fid,*,iostat=i) k
   if(i /= 0) exit
   do i = 1, k+1, 1
    read(fid,'(A)') buf
   end do ! for i
   nframe = nframe + 1
   if(nframe == ith-1) exit
  end do ! for while
 end if

 read(fid,'(A)') buf
 read(fid,'(A)') buf
 do i = 1, natom, 1
  read(fid,*) elem(i), coor(1:3,i)
 end do ! for i
 close(fid)
end subroutine read_iframe_from_xyz

! read all frames from a given .xyz file
subroutine read_all_frames_from_xyz(xyzname, nframe, natom, elem, coor)
 implicit none
 integer :: i, j, fid
 integer, intent(in) :: nframe, natom
 real(kind=8), dimension(3,natom,nframe), intent(out) :: coor
 character(len=240) :: buf
 character*2, dimension(natom,nframe), intent(out) :: elem
 character*240, intent(in) :: xyzname

 coor = 0d0; elem = ' '
 open(newunit=fid,file=TRIM(xyzname),status='old',position='rewind')

 do i = 1, nframe, 1
  read(fid,'(A)') buf
  read(fid,'(A)') buf
  do j = 1, natom, 1
   read(fid,*) elem(j,i), coor(1:3,j,i)
  end do ! for j
 end do ! for i
 close(fid)
end subroutine read_all_frames_from_xyz

! read the number of frames from pdb file
subroutine read_nframe_from_pdb(pdbname, nframe)
 implicit none
 integer :: i, fid
 integer, intent(out) :: nframe
 character(len=240) :: buf
 character*240, intent(in) :: pdbname

 nframe = 1 ! initialization
 open(newunit=fid,file=TRIM(pdbname),status='old',position='append')
 do while(.true.)
  BACKSPACE(fid,iostat=i)
  if(i /= 0) exit
  BACKSPACE(fid,iostat=i)
  if(i /= 0) exit
  read(fid,'(A)') buf
  if(buf(1:5) == 'MODEL') exit
 end do ! for while

 close(fid)
 if(i /= 0) return ! assume 1 frame

 i = INDEX(buf, ' ')
 read(buf(i+1:),*) nframe
end subroutine read_nframe_from_pdb

! read the i-th frame from a given .pdb file
! Return cell, elem, resname and coor.
! If the cell size is not recorded in the pdf file, the array cell will be 0.
! If there is no residue name recorded, the array resname will be ' '
subroutine read_iframe_from_pdb(pdbname, iframe, natom, cell, elem, resname, coor)
 implicit none
 integer :: i, j, fid, iatom
 integer, intent(in) :: iframe, natom
 real(kind=8), dimension(6), intent(out) :: cell
 real(kind=8), dimension(3,natom), intent(out) :: coor
 character(len=6) :: str
 character(len=240) :: buf
 character*240, intent(in) :: pdbname
 character*2, dimension(natom), intent(out) :: elem
 character*3, dimension(natom), intent(out) :: resname

 elem = '  '    ! initialization
 resname = '   '
 cell = 0d0 ! a,b,c,alpha,beta,gama
 coor = 0d0

 if(natom <= 0) then
  write(6,'(A)') 'ERROR in subroutine read_iframe_from_pdb: natom<=0.'
  write(6,'(A,I0)') 'Your input natom=', natom
  stop
 end if

 open(newunit=fid,file=TRIM(pdbname),status='old',position='rewind')
 do while(.true.)
  read(fid,'(A)',iostat=i) buf
  if(i /= 0) exit
  if(buf(1:5) == 'MODEL') then
   read(buf(6:),*) j
   if(j == iframe) exit
  end if
 end do ! for while

 if(i /= 0) then
  if(iframe /= 1) then
   write(6,'(A)') 'ERROR in subroutine read_iframe_from_pdb: fail to read&
                  & the i-th frame in file '//TRIM(pdbname)
   write(6,'(A,I0)') 'iframe=', iframe
   close(fid)
   stop
  else ! iframe == 1
   rewind(fid)
   do while(.true.)
    read(fid,'(A)',iostat=i) buf
    if(i /= 0) exit
    if(buf(1:4)=='ATOM' .or. buf(1:6)=='HETATM') exit
   end do ! for while
   if(i /= 0) then
    write(6,'(A)') 'ERROR in subroutine read_iframe_from_pdb: failed to read&
                   & the 1st frame in file '//TRIM(pdbname)
    close(fid)
    stop
   end if
   BACKSPACE(fid)
  end if
 end if

 do i = 1, 2
  BACKSPACE(fid,iostat=j)
  if(j /= 0) exit
  BACKSPACE(fid,iostat=j)
  if(j /= 0) exit
  read(fid,'(A)') buf
  if(buf(1:6) == 'CRYST1') then
   read(buf,*) str, cell(1:6)
   exit
  end if
 end do ! for i

 do i = 1, 3
  read(fid,'(A)') buf
  if(buf(1:4)=='ATOM' .or. buf(1:6)=='HETATM') exit
 end do ! for j
 BACKSPACE(fid)

 read(fid,*) resname(1:4)
 resname(4) = ADJUSTL(resname(4))
 read(resname(4),fmt=*,iostat=i) j
 BACKSPACE(fid)
 iatom = 0

 if(i /= 0) then ! the 4-th column is residue name
  do while(.true.)
   read(fid,'(A)') buf
   if(buf(1:3) == 'TER') cycle
   iatom = iatom + 1
   read(buf,*) str, j, elem(iatom), resname(iatom), j, coor(1:3,iatom)
   if(iatom == natom) exit
  end do ! for while
 else            ! the 4-th column is integer
  do while(.true.)
   read(fid,'(A)') buf
   if(buf(1:3) == 'TER') cycle
   iatom = iatom + 1
   read(buf,*) str, j, elem(iatom), j, coor(1:3,iatom)
   if(iatom == natom) exit
  end do ! for while
 end if

 close(fid)
 do i = 1, natom, 1
  j = IACHAR(elem(i)(2:2))
  if(j>47 .and. j<58) elem(i)(2:2) = ' '
 end do ! for i
end subroutine read_iframe_from_pdb

! read Cartesian coordinates from a (Open)Molcas output file
subroutine read_coor_from_molcas_out(outname, natom, coor)
 implicit none
 integer :: i, j, k, fid
 integer, intent(in) :: natom
 real(kind=8), intent(out) :: coor(3,natom)
 character(len=6) :: str
 character(len=240) :: buf
 character(len=240), intent(in) :: outname

 buf = ' '
 coor = 0d0
 open(newunit=fid,file=TRIM(outname),status='old',position='append')

 do while(.true.)
  BACKSPACE(fid)
  BACKSPACE(fid)
  read(fid,'(A)') buf
  if(buf(4:21) == 'This run of MOLCAS') then
   write(6,'(A)') "ERROR in subroutine read_coor_from_molcas_out: failed to&
                    & find 'Cartesian coordinates in A'"
   write(6,'(A)') 'keywords in file '//TRIM(outname)
   close(fid)
   stop
  end if
  if(buf(7:32) == 'Cartesian coordinates in A') exit
 end do ! for while

 do i = 1, 3
  read(fid,'(A)') buf
 end do

 do i = 1, natom, 1
  read(fid,*,iostat=j) k, str, coor(1:3,i)
  if(j /= 0) then
   write(6,'(A)') 'ERROR in subroutine read_coor_from_molcas_out: insufficient&
                    & number of atoms in file '//TRIM(outname)
   write(6,'(2(A,I0))') 'natom=', natom, ', but broken at i=', i
   close(fid)
   stop
  end if
 end do ! for i

 close(fid)
end subroutine read_coor_from_molcas_out

! read Cartesian coordinates from a Molpro output file
subroutine read_coor_from_molpro_out(outname, natom, coor)
 implicit none
 integer :: i, fid
 integer, intent(in) :: natom
 real(kind=8), intent(out) :: coor(3,natom)
 character(len=2) :: elem = ' '
 character(len=240) :: buf
 character(len=240), intent(in) :: outname

 coor = 0d0
 open(newunit=fid,file=TRIM(outname),status='old',position='append')

 do while(.true.)
  BACKSPACE(fid)
  BACKSPACE(fid)
  read(fid,'(A)') buf
  if(buf(2:13) == 'Current geom') exit

  if(buf(2:13) == 'Primary work') then
   write(6,'(A)') "ERROR in subroutine read_coor_from_molpro_out: no '&
                   &Current geom' found in file "//TRIM(outname)
   close(fid)
   stop
  end if
 end do ! for while

 do i = 1, 3
  read(fid,'(A)') buf
 end do

 do i = 1, natom, 1
  read(fid,*) elem, coor(1:3,i)
 end do ! for i

 close(fid)
end subroutine read_coor_from_molpro_out

! read Cartesian coordinates from an ORCA .engrad file
subroutine read_coor_from_engrad(engrad, natom, coor)
 use phys_cons, only: Bohr_const
 implicit none
 integer :: i, j, fid
 integer, intent(in) :: natom
 real(kind=8), intent(out) :: coor(3,natom)
 character(len=240) :: buf
 character(len=240), intent(in) :: engrad

 coor = 0d0
 open(newunit=fid,file=TRIM(engrad),status='old',position='rewind')

 do while(.true.)
  read(fid,'(A)',iostat=i) buf
  if(i /= 0) exit
  if(buf(1:14) == '# The atomic n') exit
 end do ! for while

 read(fid,'(A)') buf
 do i = 1, natom, 1
  read(fid,*) j, coor(:,i) ! unit in Bohr
 end do ! for i

 close(fid)
 coor = coor*Bohr_const
end subroutine read_coor_from_engrad

! write/create a Gaussian .EOu file
subroutine write_EOu(EOu, e, natom, grad)
 implicit none
 integer :: fid
 integer, intent(in) :: natom
 real(kind=8), intent(in) :: e, grad(3*natom)
 character(len=720), intent(in) :: EOu

 open(newunit=fid,file=TRIM(EOu),status='replace',position='rewind')
 write(fid,'(4D20.12)') e, 0d0,0d0,0d0
 write(fid,'(3D20.12)') grad
 close(fid)
end subroutine write_EOu

! write/create a .gjf file
subroutine write_gjf(gjfname, charge, mult, natom, elem, coor)
 implicit none
 integer :: i, fid
 integer, intent(in) :: charge, mult, natom
 real(kind=8), intent(in) :: coor(3,natom)
 character(len=240) :: chkname
 character(len=2), intent(in) :: elem(natom)
 character(len=240), intent(in) :: gjfname

 call find_specified_suffix(gjfname, '.gjf', i)
 chkname = gjfname(1:i-1)//'.chk'

 open(newunit=fid,file=TRIM(gjfname),status='replace')
 write(fid,'(A)') '%chk='//TRIM(chkname)
 write(fid,'(A)') '%nprocshared=4'
 write(fid,'(A)') '%mem=2GB'
 write(fid,'(A,//,A,//,I0,1X,I0)') '#p B3LYP/6-31G(d,p) em=GD3BJ nosymm', &
                                   'Title', charge, mult
 do i = 1, natom, 1
  write(fid,'(A2,3(1X,F18.8))') elem(i), coor(:,i)
 end do ! for i

 write(fid,'(/)',advance='no')
 close(fid)
end subroutine write_gjf

! write/create a .xyz file
subroutine write_xyzfile(natom, elem, coor, xyzname)
 implicit none
 integer :: i, fid
 integer, intent(in) :: natom
!f2py intent(in) natom
 real(kind=8), dimension(3,natom), intent(in) :: coor
!f2py intent(in) coor
!f2py depend(natom) coor
 character*2, dimension(natom), intent(in) :: elem
!f2py intent(in) elem
!f2py depend(natom) elem
 character*240, intent(in) :: xyzname
!f2py intent(in) xyzname

 open(newunit=fid,file=TRIM(xyzname),status='replace')
 write(fid,'(I0)') natom
 write(fid,'(A)') 'xyz format file produced by rwgeom of MOKIT'

 do i = 1, natom, 1
  write(fid,'(A2,3(1X,F18.8))') elem(i), coor(1:3,i)
 end do ! for i

 close(fid)
end subroutine write_xyzfile

! write a frame of molecule into a given .pdb file
subroutine write_frame_into_pdb(pdbname, iframe, natom, cell, elem, resname, &
                                coor, append)
 implicit none
 integer :: i, fid
 integer, intent(in) :: iframe, natom
!f2py intent(in) iframe, natom
 character*240, intent(in) :: pdbname
!f2py intent(in) pdbname
 character*2, dimension(natom), intent(in) :: elem
!f2py intent(in) elem
!f2py depend(natom) elem
 character*3, dimension(natom), intent(in) :: resname
!f2py intent(in) resname
!f2py depend(natom) resname
 real(kind=8), intent(in) :: cell(6), coor(3,natom)
!f2py intent(in) cell, coor
!f2py depend(natom) cell, coor
 logical, intent(in) :: append
!f2py intent(in) append

 if(append) then
  open(newunit=fid,file=TRIM(pdbname),status='old',position='append')
 else
  open(newunit=fid,file=TRIM(pdbname),status='replace')
 end if

 write(fid,'(A)') 'REMARK   1 File created by rwgeom of MOKIT'
 if(ANY(cell > 1d-4)) then
  write(fid,'(A,3(1X,F8.3),3(1X,F6.2),A)') 'CRYST1',cell(1:3),cell(4:6),' P 1           1'
 end if
 if(iframe > 0) write(fid,'(A,1X,I8)') 'MODEL',iframe

 do i = 1, natom, 1
  if(LEN_TRIM(resname(i)) == 0) then
   write(fid,'(A6,I5,2X,A2,10X,I1,4X,3F8.3,A)') 'HETATM', i, elem(i), 0, &
    coor(1:3,i),'  1.00  0.00'
  else
   write(fid,'(A4,I7,2X,A2,2X,A3,5X,I1,4X,3F8.3,A)') 'ATOM', i, elem(i), &
    resname(i), 0, coor(1:3,i), '  1.00  0.00'
  end if
 end do ! for i

 write(fid,'(A)') 'END'
 close(fid)
end subroutine write_frame_into_pdb

! calculate an internal coordinate (bond, angle, or dihedral)
function calc_an_int_coor(n, coor) result(val)
 implicit none
 integer, intent(in) :: n
 real(kind=8) :: val, rtmp(3)
 real(kind=8), intent(in) :: coor(3,n)

 val = 0d0
 select case(n)
 case(2) ! bond
  rtmp = coor(:,1) - coor(:,2)
  val = DSQRT(DOT_PRODUCT(rtmp,rtmp))
 case(3) ! angle

 case(4) ! dihedral

 case default
  write(6,'(A,I0)') 'ERROR in function calc_an_int_coor: invalid n=',n
  stop
 end select
end function calc_an_int_coor

subroutine geom_h_ring(numh, d_h_h, gjfname)
 implicit none
 integer, intent(in) :: numh
!f2py intent(in) :: numh
 real(kind=8), intent(in) :: d_h_h
!f2py intent(in) :: d_h_h
 character(len=240), intent(in) :: gjfname
!f2py intent(in) :: gjfname

 write(6,'(/,A)') 'ERROR in module rwgeom: the function geom_h_ring has been re&
                  &named to gen_h_ring'
 write(6,'(A)') 'since MOKIT-1.2.6rc19. Please use gen_h_ring now.'
 stop
end subroutine geom_h_ring

! generate a Hydrogen ring
subroutine gen_h_ring(numh, d_h_h, gjfname)
 implicit none
 integer :: i, mult
 integer, intent(in) :: numh
!f2py intent(in) :: numh
 real(kind=8) :: r, theta, c
 real(kind=8), intent(in) :: d_h_h ! in Angstrom
!f2py intent(in) :: d_h_h
 real(kind=8), parameter :: PI = 4d0*DATAN(1d0)
 real(kind=8), allocatable :: coor(:,:)
 character(len=2), allocatable :: elem(:)
 character(len=240), intent(in) :: gjfname
!f2py intent(in) :: gjfname

 mult = 1
 if(MOD(numh,2) == 1) mult = 2
 allocate(coor(3,numh), source=0d0)

 theta = 2d0*PI/DBLE(numh)
 r = d_h_h*DSIN((PI-theta)*0.5d0)/DSIN(theta)

 do i = 1, numh, 1
  c = theta*DBLE(i-1)
  coor(1,i) = DCOS(c)
  coor(2,i) = DSIN(c)
 end do ! for i
 coor = coor*r

 allocate(elem(numh))
 elem = 'H '

 call write_gjf(gjfname, 0, mult, numh, elem, coor)
 deallocate(elem, coor)
end subroutine gen_h_ring

! generate a linear Hydrogen chain
subroutine gen_h_chain(numh, d_h_h, gjfname)
 implicit none
 integer :: i, mult
 integer, intent(in) :: numh
!f2py intent(in) :: numh
 real(kind=8), intent(in) :: d_h_h ! in Angstrom
!f2py intent(in) :: d_h_h
 real(kind=8), allocatable :: coor(:,:)
 character(len=2), allocatable :: elem(:)
 character(len=240), intent(in) :: gjfname
!f2py intent(in) :: gjfname

 mult = 1
 if(MOD(numh,2) == 1) mult = 2
 allocate(coor(3,numh), source=0d0)

 do i = 1, numh, 1
  coor(3,i) = DBLE(i-1)*d_h_h
 end do ! for i

 allocate(elem(numh))
 elem = 'H '

 call write_gjf(gjfname, 0, mult, numh, elem, coor)
 deallocate(elem, coor)
end subroutine gen_h_chain

! replace Cartesian coordinates in .fch(k) file by coordinates from .gjf
subroutine replace_coor_in_fch_by_gjf(gjfname, fchname)
 implicit none
 integer :: natom
 real(kind=8), allocatable :: coor(:,:)
 character(len=240), intent(in) :: gjfname, fchname
!f2py intent(in) :: gjfname, fchname

 call read_natom_from_gjf(gjfname, natom)
 allocate(coor(3,natom))
 call read_coor_from_gjf_or_xyz(gjfname, natom, coor)
 call replace_coor_in_fch(fchname, natom, coor)
 deallocate(coor)
end subroutine replace_coor_in_fch_by_gjf

! replace Cartesian coordinates in .fch(k) file by coordinates from .engrad
subroutine replace_coor_in_fch_by_engrad(engrad, fchname)
 implicit none
 integer :: natom
 real(kind=8), allocatable :: coor(:,:)
 character(len=240), intent(in) :: engrad, fchname
!f2py intent(in) :: engrad, fchname

 call read_natom_from_engrad(engrad, natom)
 allocate(coor(3,natom))
 call read_coor_from_engrad(engrad, natom, coor)
 call replace_coor_in_fch(fchname, natom, coor)
 deallocate(coor)
end subroutine replace_coor_in_fch_by_engrad

! convert a .fch(k) file into a .xyz file
subroutine fch2xyz(fchname)
 implicit none
 integer :: i, natom, charge, mult
 integer, allocatable :: nuc(:)
 real(kind=8), allocatable :: coor(:,:)
 character(len=2), allocatable :: elem(:)
 character(len=240) :: xyzname
 character(len=240), intent(in) :: fchname
!f2py intent(in) :: fchname

 call find_specified_suffix(fchname, '.fch', i)
 xyzname = fchname(1:i-1)//'.xyz'

 call read_natom_from_fch(fchname, natom)
 allocate(nuc(natom), coor(3,natom), elem(natom))
 call read_elem_and_coor_from_fch(fchname, natom, elem, nuc, coor, charge, mult)
 deallocate(nuc)

 call write_xyzfile(natom, elem, coor, xyzname)
 deallocate(elem, coor)
end subroutine fch2xyz

! convert a .gjf file into a .xyz file
subroutine gjf2xyz(gjfname)
 implicit none
 integer :: i, natom, charge, mult
 integer, allocatable :: nuc(:)
 real(kind=8), allocatable :: coor(:,:)
 character(len=2), allocatable :: elem(:)
 character(len=240) :: xyzname
 character(len=240), intent(in) :: gjfname
!f2py intent(in) :: gjfname

 call find_specified_suffix(gjfname, '.gjf', i)
 xyzname = gjfname(1:i-1)//'.xyz'

 call read_natom_from_gjf(gjfname, natom)
 allocate(elem(natom), nuc(natom), coor(3,natom))
 call read_elem_and_coor_from_gjf(gjfname, natom, elem, nuc, coor, charge, mult)
 deallocate(nuc)
 call write_xyzfile(natom, elem, coor, xyzname)
 deallocate(elem, coor)
end subroutine gjf2xyz

! convert a .xyz file into a .gjf file
subroutine xyz2gjf(xyzname)
 implicit none
 integer :: i, natom, charge, mult
 real(kind=8), allocatable :: coor(:,:)
 character(len=2), allocatable :: elem(:)
 character(len=240) :: gjfname
 character(len=240), intent(in) :: xyzname
!f2py intent(in) :: xyzname

 call find_specified_suffix(xyzname, '.xyz', i)
 gjfname = xyzname(1:i-1)//'.gjf'
 charge = 0; mult = 1  ! initialization

 call read_natom_from_xyz(xyzname, natom)
 allocate(elem(natom), coor(3,natom))
 call read_elem_and_coor_from_xyz(xyzname, natom, elem, coor)
 call write_gjf(gjfname, charge, mult, natom, elem, coor)
 deallocate(elem, coor)
end subroutine xyz2gjf

subroutine engrad2EOu(engrad, EOu)
 implicit none
 integer :: natom
 real(kind=8) :: e
 real(kind=8), allocatable :: grad(:)
 character(len=240), intent(in) :: engrad
 character(len=720), intent(in) :: EOu

 call read_natom_from_engrad(engrad, natom)
 allocate(grad(3*natom))
 call read_grad_from_engrad(engrad, natom, e, grad)
 call write_EOu(EOu, e, natom, grad)
 deallocate(grad)
end subroutine engrad2EOu

